Role of Self-Interaction Effects in the Geometry Optimization of Small Metal Clusters 
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By combining the Self-Interaction Correction (SIC) 
with pseudopotential perturbation theory, the role of self- 
interaction errors inherent to the Local Density Approxima- 
tion (LDA) to Density Functional Theory is estimated in the 
determination of ground state and low energy isomeric struc- 
tures of small metallic clusters. Its application to neutral 
sodium clusters with 8 and 20 atoms shows that the SIC pro- 
vides sizeable effects in Nag, leading to a different ordering 
of the low lying isomeric states compared with ab-initio LDA 
predictions, whereas for Na2o, the SIC effects are less pro- 
nounced, such that a quantitative agreement is achieved be- 
tween the present method and ab-initio LDA calculations. 
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The use of Density Functional Theory (DFT) , together 
with the Local (Spin) Density Approximation (L(S)DA) 
to perform structural optimization of small finite sys- 
tems as well as to predict the stability of new, finite and 
infinite, materials, has increased dramatically in recent 
years, with a widespread use in Physics and Chemistry 
0. In this framework, the problem reduces to the re- 
peated determination of the self-consistent solution of the 
Kohn-Sham equations. Since in many materials (such 
as metals and the metallic clusters considered in this 
work) the interaction between the valence electrons and 
the ionic cores is weak - due to the effective screening 
of the nuclear charge by the tightly bound core elec- 
trons - the resulting electron-ion interaction can be for- 
mulated in terms of pseudopotentials. Several ab-initio 
norm-conserving pseudopotentials are available which en- 
sure a high degree of transferability, and are obtainable 
either via analytic formulae, or their numerical tabula- 
tion requires a few CPU seconds Indeed, ab-initio 
LSDA-DFT methods, combined with norm conserving 
pseudopotentials, have been recently used to perform 
structural minimization of clusters with as many as 147 
atoms 0. 

Other ab-initio methods, such as the Self- Consistent- 
Field Configuration Interaction (SCF-CI) which typ- 
ically treat as active all the electrons of each atom, are 
more limited in cluster size, due to their higher compu- 
tational cost, and also due to basis-set and size- 
consistency [|j problems which preclude, e.g., the study 
of cluster trends as a function of size. However, these 
methods have no self- interaction problems, such as those 
well-kown in LDA-DFT pi and, for some selected clus- 



ters (e.g., Nas), lead to ground-state structures which are 
distinct from the LDA-DFT results. The energy differ- 
ences between close lying isomers of an alkali cluster are 
often associated with charge redistribution and charge 
polarization phenomena Q], which are relevant, not only 
to distinguish between isomers of the same cluster, but 
also to understand the difference between the ground- 
state geometry of isoelectronic clusters of different alkali 
species Q|. Since these are aspects in which the self- 
interaction effects play a role, it is important to have 
an estimate of the size of these errors. In this Letter this 
problem is addressed by combining two well tested meth- 
ods: Pseudopotential Perturbation Theory (PSP- 
PT) and the Self-Interaction Correction (SIC) as pro- 
posed in ref. , together with its extension to linear re- 
sponse theories, developed in ref. It will be applied 
to optimize the geometries of two small neutral sodium 
clusters, with 8 and 20 atoms, starting from geometries 
associated with different isomers of these clusters, which 
have been identified in previous LDA-DFT ab-initio cal- 
culations jl^ (for the octamer, also in SCF-CI structural 
minimization |lj]). It will be concluded that, for Na2o, 
SIC-PSP-PT provides answers in quantitative agreement 
with the results of ref. ||l2| , both for the bond- lengths and 
for the geometries of the low energy isomers. This seems 
to indicate that SIC does not play a sizeable role in the 
structural minimization of clusters with as many as 20 
sodium atoms. However, for iVag, the results show an 
energy sequence for the three lowest isomers which does 
not follow the LDA-DFT ab-initio results of ref. As 
a result, it is observed that the structure is favoured 
with respect to the D^d isomer, in agreement with the 
SCF-CI calculations of ref. |0. This, in turn, seems to 
indicate that self-interaction effects play a sizeable role 
for the smaller clusters. 

Pseudopotential perturbation theory has been applied 
to metallic clusters in ref. using LDA-DFT. Being 
an approximation to the more sophisticated (and more 
time-consuming) ab-initio methods, it has been shown 
pO| to provide answers in quantitative agreement with 
these, in particular with state-of-the-art Car-Parrinello 
structural minimization |12| . 

The SIC prescription ^ though lacking a first prin- 
ciples justification, has been applied with great success 
in different areas of physics and chemistry, including the 
study of both ground-state and excited state properties 
of small metallic clusters At present, there is no 
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first-principles workable scheme which can remedy com- 
pletely the self-interaction errors in LDA-DFT. There- 
fore, we resort to SIC in order to estimate these effects. 
In this sense, the present study, more than conclusive, 
can be of guidance for future developments. 

Because the LDA-PSP-PT results are available for 
these clusters the role of self- interaction errors can 
be estimated unambiguosly, the only possible source of 
inaccuracy being related to the SIC itself. We would like 
to point out that SIC has been used before in this context, 
but only to test its ability to predict the bond-length in 
alkali dimers In such calculations, the role of SIC 

in the chemical bonding of atoms in a cluster cannot be 
investigated. 

On the foregoing, we briefly summarize the SIC-PSP- 
PT formulation used, before discussing the results and 
writing down the conclusions. A more detailed and ex- 
tended analysis will be published elsewhere [ p^ . 

By replacing the nuclei and core electrons by pseu- 
dopotentials, leaving as active quantal particles the va- 
lence electrons, one can write, for the ionic contribution 
to the total potential, 

N 

«ex(r,i?) = ^i;p.(f-i?,) , R^{R^} . (1) 
1=1 

In this equation, N is the number of ions (also the num- 
ber of valence electrons for the neutral clusters considered 
here), and R represents a set of given ionic positions, 
which ultimately should be determined by minimizing 
the total energy of the system. Vps{f ~ Ri) represents 
the pseudopotential at point r, centered at atomic site 
Ri. We shall adopt the local pseudopotential used in 
ref. [ p^ , with the parametrization adequate for sodium, 
taken from ref. With respect to any chosen center 

of the cluster, the sum of the pseudopotentials can be 
expanded as follows: 

Vex{r,R)= v{f,R) + V2,ex{r,R) 

= voir, R)YoAr) 

OO / 

+ ^i,m{r,R)Yi,-m{.r) . (2) 

l—l m— — l 

The first term is just the monopole, spherical part of 
the total ionic contribution, which will be taken into 
account exactly in SIC-LDA by solving the set of SIC 
Kohn-Sham-like equations for the N valence elec- 

trons moving selfconsistently in this spherical potential 
pTf . The remaining terms are included perturbatively 

up to second order, which leads to a contribution AEps 
to the total energy reading , 

AE^f - i y" drSn2{f, R) «2,ex(r, R) , (3) 

where dn2 (r, R) is the screened induced density change 
caused by the external potential V2^ex 

(f,i?). The total 



energy is then given by the sum of the monopole part, 
obtained via the SIC-LDA self-consistent solution, plus 
the term AEp'i'^ given by eq.(3). As in ref. the 
screened induced density will be computed making use 
of linear response theory, which incorporates, in a self- 
consistent manner, the screening of the external pertur- 
bation by the valence electrons. Since Sn2{r,R) is built 
out of SIC-LDA spherical self-consistent solution, it is in 
fact computed in a numerically exact way via Green's 
functions techniques (for an elaborated discussion of the 
methods and techniques, the reader is referred to refs. 

As shown in ref. the standard LDA-DFT 
linear response theory, known as Time-Dependent LDA 
(TDLDA), introduces a spurious self-polarization effect, 
which leads to an erroneous overscreening of the external 
field. The extension of SIC to the linear response theory 
has been formulated in ref. ||ll[], and it is this SIC re- 
sponse formulation which is used in computing dn2 {r, R) 

Instead of performing a fully-relaxed structural op- 
timization, we started from the given LDA-DFT low- 
energy isomeric structures found in ref. [p^ , and then 
minimize the total energy using SIC-PSP-PT by scaling 
uniformly the radial distance of the constituent atoms 
from the center of mass of the cluster. This strategy 
(similar to those adopted in refs. being incom- 

plete, will prove sufficient for the present purpose. In 
table 1 are given the results obtained with this method- 
ology for Na2o- Two low energy isomers have been iden- 
tified in ref. labelled A and B, respectively. The A 
structure corresponds to the lowest isomer, with an extra 
binding with respect to the B structure of 0.136 eV 
In the first row the PSP-PT results obtained in ref. ||lC| ] 
are shown for these 2 isomers, within LDA-DFT, using 
the minimization strategy adopted here. In the second 
row, the SIC-PSP-PT results of the present calculation 
are tabulated. All the ingredients of the 2 calculations 
are the same, except for SIC. Furthermore, and for each 
rowm the energy differences are always referring to the 
structure which has minimum energy. Furthermore, the 
numbers in parenthesis indicate the value of the scaling 
factor at which a given structure minimizes its energy. 
One can observe that SIC leads to a (small) reduction of 
the bond length in comparison with PSP-PT, in agree- 
ment with the previous results for dimers @Jl^. More 
important, however, is the fact that SIC-PSP-PT keeps 
the same energy ordering for the 2 isomers, therefore ev- 
idencing no measurable effect, apart from the expected 
difference in the restoring force coefficient for the breath- 
ing mode . 

The situation is quite different for Na^. The results 
are shown in table 2, where we keep the same notation 
as in table 1. The 3 structures considered are identi- 
fied in fig.l. Whereas SIC-PSP-PT leads to structures in 
which the volume of delocalization of the valence elec- 
trons, controlled by the scaling factor, is in excellent 
agreement with the results of ref. ||l^, the energy or- 
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dering of the isomers is different in SIC-PSP-PT, when 
compared with the predictions of refs. |]l0| , ^ . This, in 
view of the ingredients of the calculations used, turns 
out to be a pure self-interaction effect. Interestingly, the 
results obtained with the present method show a favour- 
ing of geometries which correlates well with the results 
of rcf. which makes use of a truUy self- interaction 
free theory. This, in turn, reinforces the results of the 
present calculation, which shows that self-interaction ef- 
fects should not be overlooked in the important issue of 
predicting the lowest isomeric structures of small alkali 
clusters. In this context, it is worth mentioning the re- 
cent measurements of the photoabsorption cross-section 
ofNa^ clusters Ull in a cold beam, which reveals a multi- 



peak structure which can be reproduceded |19| by com- 
puting the absorption from the ground-state geometry 
found with the self- interaction free SCF-CI after an 
appropriate scaling of the bond-lengths. 

Finally, we would like to comment on the average 
bond-length obtained using the present calculation. The 
present results lead to bond-lengths in quantitative agree- 
ment with the findings in refs. pO| , p^ , as opposed to the 
results of SCF-CI |], which lead to an average bond- 
length in small alkali clusters which is larger, on the av- 
erage, that what is known in the bulk limit As the 
cluster size decreases, one expects a contraction of the 
average bond-length, due to the corresponding increase 
of the surface to volume ratio. This qualitative reason- 
ing, as well as the comparative performance of LDA-DFT 
methods and SCF-CI methods for the sodium dimer (cf. 
ref. iQ) lend support to the results of the present cal- 
culation. In any case, the reason behind this difference 
cannot be attributed to a self-interaction problem but, in- 
stead, to the different treatment of correlations between 
the active electrons. 

In summary, in what concerns the geometry optimiza- 
tion of small alkali metal clusters, one has shown that 
self-interaction effects cannot be overlooked, and may be 
crucial for an accurate assignment and sorting of the low 
energy isomers. This feature decreases with increasing 
cluster size, and seems to be unimportant already for 
Na2o- As is well known self-interaction effects still 
manifest themselves in other properties of these clusters, 
such as the energy of the one-electron orbitals, modifying 
also the topology of the potential energy surface, which 
is reflected in a different vibrational spectrum p^ . 
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TABLE I. Results for Na2o, for the two lowest isomers 
found in ref. jl^, obtained using standard PSP-PT [|l^, and 
the SIC-PSP-PT developed in this work. For each of the 2 
geometries selected, two sets of values are given, associated 
to the result of a total energy minimization as a function of 
the scaling parameter. The value of the scaling parameter 
at minimum is given in parenthesis, whereas the value of the 
total energy constitutes the other result tabulated. The values 
quoted for the energies (eV) are referred to the state which, 
for each method, corresponds to the absolute minimum. 



Na2o 


A 


B 


PSP-PT 


0.000(1.02) 


0.111(1.01) 


SIC-PSP-PT 


0.000(1.00) 


0.237(1.00) 
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TABLE II. Results for Nag, for the three lowest isomers 
found in refs. illustrated in fig.l. The notation used 

here is the same as in table 1. 



Nas 


D2d 


Td 


Did 


PSP-PT 


0.000(1.02) 


0.159(1.02) 


0.111(1.01) 


SIC-PSP-PT 


0.000(1.00) 


0.104(1.00) 


0.159(1.00) 



Figure available from: ekardt@fhi-berlin.mpg.de 



FIG. 1. The 3 different geometries obtained in ref. |l^] 
for Nas, corresponding to the ground-state and two lowest 
isomers (energy is increasing from left to right), and used as 
starting geometries which are minimized as a function of a 
dimensionless scaling parameter, which uniformly and simul- 
taneously changes the radial distances of all the atoms with 
respect to the center of mass of the cluster. The structure 
on the left displays D2d symmetry, the one in the middle Did 
symmetry, and the one on the right Td symmetry (see main 
text for details). 
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